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Abstract 

The rosette-shaped motion of a particle in a central force field is 
known to be classically solvable by quadratures. We present a new 
approach of describing and characterizing such motion based on the 
eccentricity vector of the two body problem. In general, this vector is 
not an integral of motion. However, the orbital motion, when viewed 
from the nonuniformly rotating frame defined by the orientation of 
the eccentricity vector, can be solved analytically and will either be 
a closed periodic circulation or libration. The motion with respect to 
inertial space is then given by integrating the argument of periapsis 
with respect to time. Finally we will apply the decomposition to a 
modern central potential, the spherical Hernquist-Newton potential, 
which models dark matter halos of galaxies with central black holes. 

1 Introduction 

1.1 Central Force Fields 

The motion of a particle in a central force field is known to be classically 
solvable by quadratures. Due to the spherical symmetry of the force field, 
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an angular momentum integral exists and the ensuing motion is confined to 
a single orbital plane so that, without loss of generality, we can assume the 
system to have 2 degrees of freedom. In polar coordinates, the Hamiltonian 
may be expressed as: 

where h = is the angular momentum and U{r) is the potential energy 
function. We will use the standard convention that dots refer to time deriva- 
tives, whereas primes refer to spatial derivatives. The corresponding Hamil- 
tonian system, 

h'^ . h ■ 

r = Vr Vr = ^ + U'{r) = — h = 0, 



is integrable by quadratures: 

±dr 



rio) y/2H + 2U{r)-hyr 
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t - 1„ (1) 



9W=»o+/"'^. (2) 



to 



The ensuing motion follows rosette-shaped paths (Arnold (1989), Whittaker 
(1988), etc.). 



1.2 Osculating Orbital Elements 

One could, alternatively, proceed using Variation of Parameters and La- 
grange's Planetary Equations (Brouwer & Clemence (1961), Roy 1988). In 
this case, one can write down differential equations of motion for the six os- 
culating classical orbital elements and then solve them by quadrature for aU 
time. These equations are nonlinear and furthermore depend upon a choice 
of the "planetary" gravitational parameter /i. One assumes the potential is 
a perturbation of a Newtonian potential: 

U{r) = ^ + R{r), 
r 

where R{r) is known as the disturbing function. For a general central force 
field where there is no nominal attracting body, such a choice is somewhat ar- 
bitrary. For instance, there is no primary "planet" when considering motion 
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in a galactic halo. A gravitational parameter for the unperturbed motion can 
nonetheless be artificially contrived, perhaps based on the total halo mass (if 
the motion evolves in the outskirts of the galaxy) or based on the mass of 
a central galactic bulge or black hole (if the motion evolves near the galac- 
tic core). Whatever the choice of gravitational parameter, a complete set of 
osculating orbital elements arises and the ensuing motion can be determined. 

1.3 The Eccentric Frame 

We will define a gravitational parameter based on the central force field's 
potential, with no reference to a main attracting body and perturbation 
theory. Following the classical analogy, we define /i(r) such that: 



We will show that this gives rise to a nonstatic eccentricity vector that rotates 
at a nonuniform rate. The eccentricity vector (Runge-Lenz vector) associ- 
ated with this spatially variable gravitational parameter function defines a 
preferred coordinate system which we call the eccentric frame. With respect 
to this frame, we will show that the motion follows a closed orbit. Depending 
on the value of energy, the particle will make closed circulations or librations 
in the eccentric frame. The eccentric frame decomposition gives rise to a set 
of orbital elements. We will discuss their physical implications and the key 
features of how they arise. In particular, one can have circular orbits in iner- 
tial space with nonzero osculating eccentricity. This feature is not unique to 
our method, it can arise from any choice of osculating orbital elements. The 
eccentric frame decomposition, however, illuminates the behavior and gives 
rise to a new standard description that better fits orbits of central force field 
potentials. 

2 The Eccentric Frame Decomposition 

We first define the eccentric frame by means of specifying the nonstatic ec- 
centricity vector associated with the gravitational parameter function /i(r) = 
rU{r). We then show that the particle traces a closed orbit as viewed from 
this noninertial frame. Finally we compute the set of osculating orbital ele- 
ments that belong to this system. 




r 
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2.1 Motion with respect to the Eccentric Frame 



Given a spherically symmetric potential energy field, we can recast the Hamil- 
tonian into the following form, reminiscent of its classical analogy: 

2 \ / r 



where h = r^9 is the magnitude of the angular momentum vector, 

H = r X r, (4) 

and (r, v, 9, h) are the symplectic coordinates, with v = r. This gives rise to 
the following Hamiltonian equations of motion: 



r = V V = — 

6 = ^ h = 

which can be recast in the following form 

- ^1 e. = e,, (5) 

where r = re^. Consider now the eccentricity vector 

B = r X H — ^{r)er. (6) 
Its evolution is governed by the following equations of motion: 
B = r X H — jj!{r)rer — fi{r)9eg 



/i'(r) fx{r] 



H'{r) fi{r) 



,2 



r^9er X H — fi'(r)rer — fi{r)9eg 
r'^9eQ — fj,'{r)rer — fj,{r)9eQ 



= — /i'(r)r. 

The vector B itself works out to be 

B = r'^9{rer + r^e^) x H — /i(r)er 
= —r^r9e0 + r^9'^er — yu(r)er 

= -hree + jj{r] 

\ r 



We thus find the magnitude of B is: 



(7) 



We define the argument of periapsis, u, to be the angle made between the 
inertial x-axis and the B-vector. The B-vector defines a rotating reference 
frame, which we call the eccentric frame. We define B to be a unit vector in 
the B direction. Hats will denote unit vectors. Let X and Y be the cartesian 
coordinates of the particle with respect to the eccentric frame and let x and 
y be the cartesian coordinates of the particle with respect to the inertial 
frame. The axes of the inertial frame are determined by the stationary unit 
vectors i and j. The polar angle of the particle measured with respect to the 
B direction is known as the true anomaly /. This notation is also used in 
Roy (1988). The polar angle of the particle in the inertial frame is related to 
the true anomaly by the following relation: 



e = f + iu. 



(8) 



In some of the literature, the true anomaly / is denoted by z/; and the 
inertial polar angle (argument of latitude) 9 is denoted by u. Decomposing 
the eccentricity vector B in the inertial frame, we see that 



B 



rh sin 6 



— /i(r) J cos6' 1 i 



-rh cos 9 



h^ 

/i(r^ 

r 



sin 6 



However, by definition, B = B{r)[cosuji + sincuj). Hence: 

1 



cosu; 
sin Co" 



B{r) 
1 

B{r) 



A 
A F 



cos 6' 
sin 6 



coscu 
sinu; 



(9) 



where we have defined 
A = 

and 



[h'^/r — /i(r)) rh 

—rh {h"^ / r — fi{r)) 



cos / — sin / 
sin / cos / 
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and have further made use of the trigonometric identities 



cos6' = cos(/ + u;) = cos / coscj — sin / sinu; 
sin6' = sin(/ + cu) = sin/coscu + cos/sino;. 

We recognize that the matrix premuhiplying the vector (cosu;, sinu;) on the 
right hand side of must be the identity matrix. Hence we have found an 
exphcit expression relating the true anomaly and the radius: 

sin/ = r^rh. (11) 

Thus we see that the particle traces out a closed path in the eccentric frame. 
By carefully considering flTT]) . we see that periapsis is always achieved at 
/ = 0, i.e. when r and B are parallel; and that apoapsis is achieved at 
/ = TT, i.e. when r and B are anti-parallel. 

If the angular momentum is positive, (fTTj) tells us that r is increasing when 
the particle is in the upper half plane and is decreasing when the particle 
is in the lower half plane. The opposite is true for the case of a negative 
angular momentum. 

2.2 The Osculating Eccentricity and Semi-Major Axis 

We would also like to point out that one can rearrange flTU]) into the following 
form: 

^ _ hVl^jr) ^ pjr) 

1 + B{r) cos f / fi{r) l + e(r)cos/' 

completely analogous to its classical (/i(r) = const.) form. 

Utilizing the relation p = a(l — e^), we can define the osculating eccen- 
tricity and semi-major axis of the system in closed form as follows: 

These are given unambiguously as a function of r, without integrating. 
They represent a standard decomposition of the motion. Using a standard 
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choice of osculating orbital elements, one would first define a semi-arbitrary 
choice for a fixed /i. Thus, there is no unique standard set of osculating 
orbital elements for a general system, rather a one parameter family of orbital 
elements that describe the motion. By using the radially varying /i(r), we 
seek to better normalize the description of motion in such systems. 

Since the true anomaly is given by the relations flTUl) and flTTl) . one now 
only need solve for the osculating argument of periapsis to obtain the com- 
plete motion as a function of time. 

2.3 The Osculating Argument of Periapsis 

Solving for the osculating argument of periapsis can be done in one of two 
ways. First, one may integrate ([I])- ([2D by quadratures. Once r and 9 are 
known, /, a, and e can be extracted by the above relations ffTU]l . ffTTl) . f[T^ . 
( IT^ : then the osculating argument of periapsis can be solved by means of 
the relations 6 = f + uj. On the other hand one can solve the quadrature we 
derive below. 

To determine the rotation of the eccentric frame, consider the angular 
momentum integral: 



h = r^e = r\f + dj) 



r"^ f'{r)f j 1 + 




) 



(15) 



Differentiating (fTOi) and utilizing (fTTl) . we have that 




We can now solve (fTSll for uo' 
duj 



(/): 



h 




(16) 




- 1 



$(r) = ^{r)r'^B'{r) - B'{r)rh^ 



ji' {r)r'^ B {r) . 



(17) 
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We thus have 

where we recognize r = r{f) in the integrand, by the relations ( |TOl) and ( ITT]) . 
Together with ( |T3l) and ( fT4l) . this constitutes a full set of osculating orbital 
elements that are well-defined for the orbit for all time. 

The full motion is then completely specified in terms of the parameter / 
by the relation: 

eif) = f + iu{f). 



3 The Zero Velocity Curve 

The central force problem is a 2 degree of freedom problem with 2 integrals of 
motion, E and h. It is therefore integrable and, in fact, reduces to motion on a 
Liouville torus. The symplectic coordinates of the system are (r, v, 6, h). The 
coordinate h is conserved, and the motion therefore takes place on the h = 
const, hyper-plane. Motion in the (r, v) plane is constrained to the curve Th^E 
defined by ([3]), with fixed E and h. Meanwhile, 6 cycles along according 
to h = Motion in the reduced (r, v, 9) space can therefore be visualized 
as follows: it is constrained to the surface obtained by revolving the curve 
Th,E around the v axis. This resulting surface is (obviously) topologically 
equivalent to the Liouville torus, but obtained directly without the arduous 
task of computing action-angle variables. 

For a fixed as one varies the energy, one encounters various bifurcation 
points where the system undergoes changes. 



3.1 Periapsis and Apoapsis 

Computation of the periapsis and apoapsis radii is accomplished by the stan- 
dard technique of plotting the zero-velocity curve on the E — r plane. The 
plot is obtained by setting f = r = in ([3]). The resulting equation is: 



^.vM4^-^. (19) 



For a fixed energy E, the solutions to this equation represent the periapsis Vp 
and apoapsis radii. Maximum and minimum values oi oi E^^ir) correspond 
to unstable and stable circular orbits, respectively. If there are multiple 
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"wells," the corresponding roots of this equation alternate r^i, Vai, rp2,ra2, ■ ■ 
and the forbidden regions of the inertial x-y plane are concentric, circular 
annuli. 

3.2 Circular Orbits 

As one decreases the energy for a fixed angular momentum, the curves Th,E 
on the (r, ?;)-plane shrink until they degenerate to a single point on the r-axis 
which corresponds to a circular orbit in the (r, 6) polar plane. This occurs 
at the local minima of E on the (r, v) plane, and hence is given by VE = 0, 
where E is given by ([3]) and h is held fixed. This condition amounts to 

V = (20) 
ft^^^_M^ = = f. (21) 

The root of (12T]) . rcirc, corresponds to the radius of the circular orbit which 
occurs at the minimum energy i^circ := -£'zv(^circ)- 

3.3 Escape Orbits 

If U{r) const, as r ^ oo, a series of unbounded orbits are present in 
the solution space. Such orbits are classified as escape orbits. Typically one 
takes the potential at infinity to be zero, so that U{r) — > as r — *• oo, so 
that orbits with negative energies are gravitationally bounded to the center 
of the potential, whereas orbits with positive energies have enough energy to 
escape to infinity. 

4 Circulations vs. Librations in the Eccentric 
Frame 

As one decreases the energy from i^osc to the minimum energy -Ecirc, one 
encounters a bifurcation in the eccentric frame at -Ecrit, as the orbits (as seen 
from the eccentric frame) change from circulations to librations. This is a 
necessary transition that must occur, as one lowers the energy, before one 
can reach a circular orbit. It will be our goal in this section to understand 
the how this bifurcation comes about and to give a qualitative description of 
motion in the eccentric frame for fixed h as one varies E. 
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4.1 The Critical Energy 



We now define a critical radius and critical energy. The critical radius is 
defined as the root to the right hand side of ffTOl) . which occurs when: 



For a fixed h, let the solution to fl22|) be rcrit- Further, let us define the 
critical energy as follows: 



where E^^[r) is given by (fTOl) . In the following subsections, we will see how 
passing through this value of energy brings about a bifurcation in our system. 

The critical radius rcrit has an important physical signifigance in terms 
of the eccentric frame. From ( fTOj) . we see that cos / = 0, i.e. the particle 
is crossing the F-axis in the eccentric frame, exactly when r = rent- It is 
interesting to note that rcrit is independent of the energy of the system. Thus, 
as one changes the energy, the particle passes between the left and right hand 
planes through the same two portals {Y = ±rcrit)- 

As one decreases the energy, the zero velocity curves r = Tp and r = r^ 
come closer together. Eventually, one will coincide with rcrit- This occurs 
at the critical energy E^nt and brings about the bifurcation in the system. 
For E < -Ecrit, the points Y = ±rcrit both lie in the forbidden region, thus 
a transition from the left half to right half plane is no longer possible. If 
the apoapsis zero-velocity curve r = r^ reaches rcrit before the periapsis 
zero-velocity curve r = Vp does, the particle follows periapsis librations (i.e. 
librations around periapsis on the right half plane) in the eccentric frame. 
Alternatively, if the periapsis zero- velocity curve r^ reaches rcrit first, the 
particle follows apoapsis librations. 

4.2 The Route to Periapsis Librations 

We will first consider the case where it is the apoapsis radius that coincides 
with rcrit at the bifurcation energy -Ecrit- This event brings about periapsis 
librations for all energies E < -Ecrit- In Fig- [U and Fig. [2l the path of the 
particle, for various values of energy, is plotted with respect to the eccentric 
frame, i.e. the X-axis is coincident with the eccentricity vector B. As viewed 




(22) 




1 At(^crit) 



(23) 
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from this nonuniformly rotating frame, the trajectory of the particle makes 
closed orbits. 

If £■ >> -Ecrit, fflOl) and ffTTl) produce a well-defined closed orbit in the 
eccentric frame, as seen in Figured^. 

As E approaches E'cnt from above, the apoapsis radius slowly approaches 
the critical radius, and an orbit such as the one seen in Figure [T]d is present. 
Notice the left half of this orbit is nearly circular. This presents some nu- 
merical difficulty if one discretizes the radius r and not the true anomaly 
/. However, this difficulty can be overcome by analytically approximating 
the left half of the orbit with an ellipse rapprox(/), / £ [~7r/2, — vr] U [vr/2, ir] 
fitted to the data points rapprox(±7r/2) = r^it and rapprox(vr) = r^. 

At E = Ecrit, the apoapsis radius and the critical radius coincide, as 
shown in Figure [T]3. The particle thus reaches the F-axis of the eccentric 
frame at the precise moment it reaches the zero velocity curve. Recall that 
ffTTj) implies that r is increasing in the upper half plane and decreasing in 
the lower half plane for the case h > 0. This bifurcation point is rather 
interesting, as one only has a half orbit in the eccentric frame. The motion 
begins at periapsis, but when it reaches the F-axis, i.e. apoapsis, it "hops" vr- 
radians to the corresponding point on the lower half plane and then returns 
to periapsis. To compensate there is a corresponding vr-radian hop in the 
argument of periapsis, so that the true polar angle 6' is a continuous function 
of time. This is allowed as the ^ equation, (1161) . is actually undefined for 
r = 0. This is permissible because B{rcrit) = exactly if E = Ecnt, i-e. the 
eccentricity vector actually vanishes at these endpoints, and then reappears 
pointing in the opposite direction. 

For E < Ecnt, one sees that < rcrit- Thus the portal X = 0, Y = ±rcrit 
from the right- to the left- half plane lies in the forbidden region. Motion is 
therefore constrained to the right-half, where periapsis librations arise in the 
eccentric frame. Figure [T]i-[I^. 

Finally, &t E = Edrc, the periapsis and apoapsis radii coincide and the 
trajectory in the eccentric frame degenerates to a single point X = r^nt, Y = 
0, as seen in Figure [1}:. The eccentric frame now rotates at a uniform rate 
and a circular orbit is present in the actual inertial space. 

4.3 The Route to Apoapsis Librations 

As one sees from Figure [2l the case where the periapsis radius and the critical 
radius coincide at the bifurcation energy E = E^nt leads to apoapsis librations 
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Figure 1: The Route to Periapsis Librations 



in the left half plane. 

An easy test to determine whether the librations will be periapsis or 
apoapsis librations is as follows: 

'"circ < '''crit =^ periapsis librations 
Tcirc > ^crit =^ apoapsis librations. 

Again one sees that there is a vr-radian hop in both true anomaly / and 
argument of periapsis u at the bifurcation energy E = E^nt- Just before the 
bifurcation, the right half of the orbit (the half closest to periapsis) is nearly 
circular. 



5 Symmetry of the Rotation 

One can exploit the form of the dynamical equation for u (1161) to reduce the 
numerical integration to one over only one half of an orbit. By examining the 
differential equation f|T6l) . one sees that uj'{f) depends only upon the radial 
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Figure 2: The Route to Apoapsis Librations 

coordinate r. Due to the periodicity of the orbit, we have that uj'{f) is 27r 
periodic. Moreover, for / G [7r,27r], we have that uj'{f) = uj{27i — /), since 
the orbits in the eccentric frame are symmetric with respect to the x-axis. 

5.1 Circulations 

Prior to the bifurcation [E > -Ecrit) the trajectory makes closed circulations 
in the eccentric frame. During the circulations, there is a secular growth in 
the argument of periapsis u. Define the following: 

n = ^{r{f)) df /G[0,7r], 

such that Q is the argument of periapsis restricted to the domain / G [0, vr]. 
We will show that once one has fl, one can systematically find uj{f) for all 
future /, without integration. 

The condition that uj'{f) = uj'{2'K — f) for / G [vr, 27r] suggests that the 
function uj{f) is odd with respect to the axes / = tt and u = uj{tt) on 
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the interval [0,27r]. Thus, given Cl (which we presume has been found by a 
numerical algorithm), one defines 

0(f) = ! ^(/) /e[0,7r] 

\ 2n{n)-n{2n-f) / G (tt, 27r] " 

The net secular growth in uj{f) over one nominal orbit < / < 27r is given 
by 

An = 0(27r). 

uj{f) can subsequently be found by applying the following: 

uj{f) = nAQ + n{f mod 27r), 

where n is the orbit number, i.e. n = if / G [0, 2n], n = 1 if / G [27r, 47r], 
etc. 

6 The Hernquist-Newton Potential 

To illustrate the theory in the context of a modern problem, we will consider 
motion of a particle (star) in a spherical galaxy, modelled with the Hern- 
quist potential, with a central black hole. These results could be similarly 
applied to a black hole at the center of a globular cluster, or various other 
astrophysical configurations that yield spherical or azimuthal symmetry. In 
this context, the central black hole provides a classical point potential, but 
no general relativistic effects are included. 

6.1 Galactic Halos with Central Black Holes 

The Hernquist potential has achieved some acclaim in recent years for its 
ability to analytically model galactic dark matter halos, see Hernquist (1990). 
We will consider here a coupling between the spherical Hernquist profile and 
a Newtonian point mass, assumed to model a black hole at the center of 
the galaxy. Some numerical modelling of triaxial galaxies with central black 
holes has already been carried out, as in Poon & Merrit (2004). 

Let /iBH and yUhaio be the gravitational parameters of the central black 
hole and the galactic dark matter halo, respectively; and let 6 be a length 
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scale of the galaxy (so that M{b) = Mtot/4, see Hernquist (1990)). Then the 
Hernquist-Newton potential can be written: 



U{r) 



/^halo /^BH 

R + b ~W 



By defining: 



fJ'O — fJ'halo + A'BH fJ' 



/^halo 



/^halo + /^BH 

the Hernquist-Newton potential can be recast into the following equivalent 
form: 

U{r) = 

with associated Hamiltonian: 



1 



R/b 



dR 

It 



IP 

i?2 



R 



R/b 



(24) 



where H = R^^ is the angular momentum. As this is a central force field, 
the angular momentum and energy will be conserved quantities. Observe 
that when fi = 0, the potential energy reduces to that of a Newtonian point 
mass. When fl = 1, the potential energy is equivalent to the Hernquist 
potential. For < /i << 1, the model represents a Newtonian point mass 
with a surrounding "Hernquist cloud" and for << /2 < 1, we have the 
Hernquist potential with a relatively weak point mass at the origin, which 
could be used to model a spherical Hernquist galaxy with a central black 
hole. 



6.2 Nondimensionalization 

Carrying out the following change of variables: 



and thus, consequently. 
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we can recast the Hamiltonian fl24p into the following form: 



where we define 



/^(^) = ( 1 - ) • (26) 



We have thus recast the Hernquist-Newton potential to a one-parameter 
family of potentials, with fl = 1 corresponding to a the Hernquist potential 
and jl = corresponding to a pure Newtonian point mass. 

We note that the Hernquist-Newton potential is similar to analogous work 
on the Manev problem, which considers a potential of the form U{r) = 
A/r + B/r"^. In fact, work has been carried out for the anisotropic Manev 
problem, which replaces the radial coordinate r with an "elliptic radius" 
m = ^y ij,x'^ + y"^ (e.g., Craig et al. 1999, Diacu & Santoprete 2001). In 
this type of potential, one obtains a large class of chaotic orbits as well as 
nonchaotic orbits. 



6.3 Zero Velocity Curves 

Using the relationship for circular orbits fl2T|) and the critical energy condition 
( 122|) . sample values of circular radius and energy, critical energy, and the 
critical periapsis and apoapsis are shown in Table [H where we have taken 
h = 0.1. The critical radius r^rit coincides with the critical energy periapsis 
radius r^°" in each case, so that the bifurcation always leads to apoapsis 
librations. 



fl 


'"circ 


E ■ 


-E'crit 


pen 
^crit 


apo 
^crit 


1 


0.2500 


-0.7200 


-0.4524 


0.1051 


1.1932 


0.99 


0.2274 


-0.7539 


-0.5000 


0.1000 


1.000 


0.95 


0.1508 


-0.9372 


-0.7440 


0.0820 


0.4460 


0.90 


0.0938 


-1.3206 


-1.1960 


0.0647 


0.1645 



Table 1: Various physical quantities for /i = 1 



The zero- velocity curves are plotted in Fig. [3] below for h = 0.1 and for 
the same values of fl as in Table [H 
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2 



r 



Figure 3: Zero Velocity Curves with h = 1 and, from top down, = 
1,0.99,0.95,0.9 

6.4 Orbits for h = 0.1 



We examine a sample of orbits in the eccentric and inertial frames for various 
energies at h = 0.1. The zero velocity curve for this angular momentum is 
plotted in Fig. |4] below. 

For a sample orbit with E > -Ecrit, we take E = —0.6. The orbit as seen 
from the eccentric and inertial frames is shown in Fig. [51 

Upon integrating fllSp . one obtains uj{f), which can be seen for this orbit 
plotted in Fig. O One sees Au = a;(27r) — ^^(O) is the turning angle of 
the rosette. For energies prior to (above) the critical bifurcation energy, we 
find a secular retrograde rotation of the eccentric frame. 0{t), uj{t), and f{t) 
are plotted against time over three standard orbits on the right. One sees 
a secular prograde growth in the inertial polar angle 6, with turning angle 
Ae = 27T + Au}. (Recall that Au <0). 

Finally, we compute the osculating semi-major axis and eccentricity vs. 
time (Fig. [7]) over one nominal orbit (as seen from the eccentric frame). 
The solid lines represent the osculating elements as provided by the eccentric 
frame method, see (fT^ - (IT^ . The dashed curves are a standard set using the 
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0.4 



h=0.1 zero velocity curve 



0.2- 
0- 




_|i , , , 1 

0.5 1 1.5 2 



r 

Figure 4: Zero Velocity Curves and sample orbits ai h = 0.1; E = 

— 0.6, -Ecrit, —0.8, -Ecirc 




Figure 5: E = —0.6 orbit in eccentric (left) and inertial (right) frames 



osculating orbital element transformation as defined by classical perturbation 
theory, using /ig = 1 for the "planet" mass, i.e. /io is the gravitational 
parameter of the total halo mass plus the mass of the central black hole. 
For a sample orbit with E < Ecnt, we take E = —0.8. The orbit as seen 
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-1.5 




0.25 

o 

■| 0.2 
io.15 
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0.05 



- polar angle; 6 

- - true anomaly: t 

- - argument of periapsis; to 



Figure 6: io{f) (left) and e{t), uj{t), f{t) (right) for E = -0.6 



various osculating semi-major axis 



various osculating eccentricity 
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Figure 7: Osculating orbital elements for E = —0.6 



from the eccentric and inertial frames is shown in Fig. [8l Notice that the 
particle now makes librations in the eccentric frame. The libration paths 
become smaller and smaller until they degenerate to a single point at Echc, 
to be considered next. 

uj{f) for this orbit is plotted in Fig. [91 Notice that / librates around 
/ = vr and there is secular prograde growth in u. The turning angle is still 
given by Au. For energies after (below) the critical bifurcation energy, we 
find a secular prograde rotation of the eccentric frame. 0(t), uj(t), and f{t) 
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Figure 8: E = —0.8 orbit in eccentric (left) and inertial (right) frames 



are plotted against time over three standard orbits on the right. One now 
sees a secular posigrade growth in the argument of periapsis, coupled with 
librations in true anomaly /. 




f t 
Figure 9: uj{f) (left) and e{t), uj{t), f{t) (right) for E = -0.8 



Finally, we compute the osculating semi-major axis and eccentricity vs. 
time (Fig. [TOl) over one nominal orbit (as seen from the eccentric frame). 
The solid lines represent the osculating elements as provided by the eccentric 
frame method, see (fT3|) - f|T^ . As before, the dashed curves are a standard set 
using the osculating orbital element transformation as defined by classical 
perturbation theory, using /iq = 1 for the "planet" mass, i.e. /iq is the 
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gravitational parameter of the total halo mass plus the mass of the central 
black hole. 




2 



Figure 10: Osculating orbital elements for E = —0.8 



At the circular energy E = Echc, the orbit in the eccentric frame degen- 
erates to a single point. The orbit is circular in the inertial plane, Fig. [TTl 
The eccentric frame now precesses at a uniform rate. 
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Figure 11: E = —0.9372 orbit in eccentric (left) and inertial (right) frames 



uj{f) (left) and e{t),u;{t)J{t) are plotted below in Fig. [H Now / is 
virtually constant and there is uniform growth in u. 

Finally, when observing the semi-major axis and eccentricity of the orbit 
(Fig. [T51) . we see something counterintuitive. The osculating eccentricity is 
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close to 0.62. If one, on the other hand, used a classical definition of oscu- 
lating orbital elements, as previously discussed, the osculating eccentricity 
would be close to 0.94. We thus see a circular orbit (in inertial space) with 
high osculating eccentricity. The osculating ellipse is a highly eccentric one 
that always touches the true path at apoapsis. In this way, the osculating 
ellipse rotates synchronously with the particle so that the particle is always 
at apoapsis and the true motion is a circular path. 
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Figure 13: Osculating orbital elements for E = —0.9372 
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7 Stability Analysis of Equatorial Rosettes in 
Axi-symmetric Potentials 



In this section we will consider the stability of the planar equatorial motion 
of a particle in an axi-symmetric potential. On the plane of symmetry, the 
potential reduces to a two-dimensional central force field. Since the motion 
is periodic with respect to the eccentric frame, the technique justifies ap- 
plication of Floquet theory for the stability analysis of the in-plane motion, 
where the period used is the period of the orbit in the eccentric frame. 



7.1 Floquet Theory 



Suppose now we are considering a conservative system with axi-symmetric 
equipotential contours, so that, in cylindrical coordinates, the potential en- 
ergy is given by U{r, z). We can define a variable gravitational parameter by 
^(r, z) = rU{r, z). The Hamiltonian of the system is thus given by: 



H 



vl + 



U{r,z) 



where the symplectic coordinates are {r,Vr,9,hz,z,Vz). The equations of 
motion are: 

hi dU(r,z) 

r = Vr Vr = ^ + 



K = o 

dU{r, z 



dr 



oz 

In particular, equatorial motion with 2; = is well defined, and this 
reduces to the problem of motion in a central force field. The resulting 
motion follows a rosette-shaped path in the equatorial plane. We now ask 
whether this motion is stable under a small out-of-plane perturbation 5z. To 
do this, we consider the 2x2 out-of-plane State Transition Matrix (STM) 
$, which, by definition, gives: 



5z{t) 
6zlt) 



Sz{0) 
6zlo) 
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The STM is determined by integrating the following differential equation: 



1 

?7.,(r,0) 



1 
1 



where the coefficient matrix is evaluated along the nominal orbit {z = 0) in 
the equatorial plane. Since the equatorial motion reduces to a central force 
field problem, there exists an eccentric frame decomposition, in which the 
motion is periodic. Let T be the period of a single orbit in the eccentric 
frame. The coefficient matrix, above, only depends on r{t), and thus it is 
T-periodic. We are therefore justified in using Floquet theory in the stability 
analysis. Let Ai,A2 be the eigenvalues of $(T). Since A1A2 = 1, either the 
eigenvalues are complex conjugates on the unit circle in the complex plane, 
or real- valued with A2 = A^^. The rosettal motion on the equatorial plane 
is therefore stable if and only if both eigenvalues are on the unit circle. A 
bifurcation from stable to unstable must occur when Ai = A2 = 1. 



7.2 Application to a toy axi-symmetric potential 

To show how this theory might be applied, we will consider the following toy 
potential: 

Here, r = x"^ + and a > is a parameter. When a < 1, the potential is 
oblate spheroidal, and when a > 1, it is prolate. This potential is motivated 
by replacing r with R in ( l26l) . but does not have direct physical significance. 
Its utility to us is only to illustrate the theory. The potential reduces to the 
Hernquist potential when restricted to the equatorial plane. The question 
now arises, for various values of the parameter a, when is the equatorial 
motion due to out of plane perturbations stable? Clearly, for a = 1, the 
motion is stable due to the angular momentum integral. To proceed, we 
compute the out of plane dynamics: 

.. _dU _ -fi{R) z fi'{R) z 
^ ~ 97 ~ R^ ^ ^ R Ra^' 
For a small perturbation 6z, we obtain: 
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The coefficient uj{t) is a function of time because we have an exphcit solution 
for r(t) for the nominal motion along the equatorial plane. The out of plane 
State Transition Matrix (STM) for {6z, 6z)'^ can then be written as: 



1 

uj{t) 



m 



1 
1 



The STM $ can now be integrated along with the nominal planar solution, 
as its dynamic coefficient matrix depends only on r(t). We can now system- 
atically integrate $ between t G [0, T], where T is the eccentric frame orbital 
period. Computing the eigenvalues of $(T) reveals the stability of the planar 
equatorial orbit. 

We followed this procedure for a sampling of different energy levels and 
axis ratios 1 : 1 : a. The result for oblate potentials is shown in Fig. [TH 
The grid points with dots correspond to parameter values of a and E for 
which the equatorial motion is unstable. For prolate potentials, the result is 
similarly depicted in Fig. [151 Notice in both plots, for a = 1, that the planar 
motion is stable for all energies. 



8 Conclusion 

In this paper we presented a preferred, nonuniformly rotating frame that 
exists for motion in any central force field, with respect to which the orbital 
motion is periodic. We showed that for high values of energy, the particle 
trajectories in the eccentric frame make circulations. However, when the 
energy drops beneath a certain critical level, the trajectories follow librations 
in the eccentric frame. This is not a true bifurcation of the system, as there 
is no distinguishable physical change when the orbits are viewed with respect 
to inertial space, yet it is a necessary transition that must occur as one nears 
the minimum circular orbit energy. For circular orbits (in inertial space), 
motion in the eccentric frame degenerates from librations to a single fixed 
point. For this case, the eccentric frame rotates at a constant rate, and a 
circular orbit in inertial space is observed. Further we showed that, even 
in the case of a circular orbit, the osculating eccentricity can be very high. 
This occurs because the particle is "stuck" to the periapsis or apoapsis of 
the osculating ellipse. The osculating ellipse has static (high) eccentricity 
and rotates at a uniform rate. We also presented a model for the potential 
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h = 0.1, |l = 1, z-stability in oblate axi-symmetric potential 
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Figure 14: Equatorial stability plot for oblate potential 

energy of a Hernquist galaxy with a central black hole, analyzed the rosette- 
shaped orbits, and then compared them to the orbit as seen from the eccentric 
frame for various parameter values. Finally we indicated how one might use 
the eccentric frame to determine stability of planar orbits in axi-symmetric 
potentials. 
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